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A-posteriori  error  estimates  for 
mixed  finite  element  and  finite  volume  methods 
for  problems  coupled  through  a  boundary 
with  non-matching  grids 

T.  Arbogast  f,  D.  EstepX ,  B.  Sheehan §  and  S.  Tavener‘S 
[Received  on  23  August  2013] 

The  primary  purpose  of  this  paper  is  to  compare  the  accuracy  and  performance  of  two  numerical  ap¬ 
proaches  to  solving  systems  of  partial  differential  equations.  These  equations  are  posed  on  adjoining 
domains  sharing  boundary  conditions  on  a  common  boundary  interface  in  the  important  case  when  the 
meshes  used  on  the  two  domains  are  non-matching  across  the  interface.  The  first  widely  used  approach 
is  based  on  a  finite  volume  method  employing  ad  hoc  projections  to  relate  approximations  on  the  two 
domains  across  the  interface.  The  second  approach  uses  the  mathematically-founded  mortar  mixed  finite 
element  method.  To  quantify  the  performance,  we  use  a  goal-oriented  a-posteriori  error  estimate  that 
quantifies  various  aspects  of  discretization  error  to  the  overall  error.  While  the  performance  difference 
may  be  not  a  surprise  in  some  cases,  we  believe  that  there  is  a  perception  in  part  of  the  scientific  commu¬ 
nity  concerned  with  multiphysics  systems  that  if  the  solution  is  smooth  near  the  interface,  then  there  is 
little  effect  from  varying  the  coupling  technique.  We  find  that,  on  the  contrary,  the  error  associated  with 
ad  hoc  coupling  approaches  may  be  large  in  practical  situations.  Moreover,  we  also  show  that  mortar 
methods  can  be  used  with  black  box  component  solves,  thus  permitting  an  efficient  and  practical  imple¬ 
mentation  also  within  legacy  codes. 

Keywords :  mortar  methods,  a-posteriori  error  estimate,  coupled  elliptic  problems,  heterogeneous  domain 
decomposition,  geometric  coupling 


1.  Introduction 

An  important  class  of  multiphysics  problems  has  a  structure  in  which  one  physical  process  dominates 
in  one  subdomain  of  the  problem  domain,  while  a  second  physical  process  dominates  in  a  neighboring 
subdomain.  The  solutions  are  coupled  by  continuity  of  state  and  continuity  of  normal  flux  through  a 
shared  boundary  between  the  subdomains.  Examples  include  general  problems  of  the  heterogeneous 
domain  decomposition  type  (Quarteroni  et  al. ,  1992;  Gaiffe  el  al.,  2002;  Bernardi  el  al.,  1994),  core- 


’  Department  of  Mathematics,  University  of  Texas  at  Austin  (arbogast@math.utexas.edu).  T.  Arbogast’s  work  was  supported  as 
part  of  the  Center  for  Frontiers  of  Subsurface  Energy  Security,  an  Energy  Frontier  Research  Center  funded  by  the  U.S.  Department 
of  Energy,  Office  of  Science,  Office  of  Basic  Energy  Sciences  under  Award  Number  DE-SC0001 1 14. 

*  Department  of  Statistics,  Colorado  State  University,  Fort  Collins,  CO  80523  (estep@stat .  colostate  .  edu).  D.  Estep’s 
work  is  supported  in  part  by  the  Defense  Threat  Reduction  Agency  (HDTRA 1-09- 1-0036),  Department  of  Energy  (DE-FG02- 
04ER25620,  DE-FG02-05ER25699,  DE-FC02-07ER54909,  DE-SC0001724,  DE-SC0005304,  INL00120133),  Idaho  National 
Laboratory  (00069249,  00115474),  Lawrence  Livermore  National  Laboratory  (B584647,  B590495),  National  Science  Founda¬ 
tion  (DMS-0107832,  DMS-0715135,  DGE-022 1595003,  MSPA-CSE-0434354,  ECCS-0700559,  DMS-1065046,  DMS-1016268, 
DMS-FRG- 1065046),  National  Institutes  of  Health  (#R01GM096192). 

§  Department  of  Mathematics,  Colorado  State  University  (brendansheehan6@gmail.com). 

^Department  of  Mathematics,  Colorado  State  University. 


(c)  The  author  2013.  Published  by  Oxford  University  Press  on  behalf  of  the  Institute  of  Mathematics  and  its  Applications.  All  rights  reserved. 


1 


2  of  24 


ARBOGAST,  ESTEP,  SHEEHAN,  TAVENER 


edge  plasma  simulations  of  a  tokamak  fusion  experiment  (Cary  et  al.,  2008,  2010),  and  conjugate  heat 
transfer  between  a  fluid  and  solid  object  (Estep  et  al.,  2008,  2009b,  2010). 

In  such  situations,  it  is  common  to  encounter  significant  differences  in  scales  of  behavior  in  the  two 
subdomains.  This  in  turn  suggests  the  use  of  different  discretization  grids.  However,  this  introduces  the 
problem  of  interpreting  the  meaning  of  coupling  state  and  flux  values  through  the  common  boundary  in 
the  discretization,  since  exact  pointwise  matching  is  no  longer  possible. 

Confounding  this  issue  are  the  practical  difficulties  of  solving  the  large  linear  and  nonlinear  discrete 
systems  associated  with  computing  numerical  solutions  and  the  common  situation  in  which  two  dif¬ 
ferent  codes  are  used  to  solve  the  two  subdomain  problems.  These  difficulties  are  generally  tackled  by 
employing  some  form  of  iterative  approach  that  involves  sequential  solution  of  the  subdomain  problems. 
The  particular  properties  of  the  discretizations  used  for  each  component  problem,  the  choice  of  iterative 
solution  method,  and  high  performance  computational  considerations  all  have  a  large  impact  on  the  way 
in  which  state  and  flux  values  are  passed  across  the  common  interface. 

In  this  paper,  we  investigate  the  accuracy  of  two  approaches  to  computing  the  coupling  values  in 
the  situation  in  which  the  discretization  grids  in  the  two  subdomains  do  not  match  at  the  interface. 
The  analysis  is  carried  out  for  the  closely  related  mixed  finite  element  and  cell-centered  finite  volume 
methods.  The  two  approaches  are  (1)  the  mortar  element  approach  (Brezzi  &  Fortin,  1991;  Roberts 
&  Thomas,  1991;  Arbogast  et  al.,  2000;  Ben  Belgacem,  2000;  Arbogast  et  al.,  2007;  Ganis  &  Yotov, 
2009),  which  uses  a  rigorous  variational  formulation  to  define  a  weak  sense  of  coupling,  and  (2)  a  “geo¬ 
metric”  approach  that  employs  various  ad  hoc  extrapolation  and  averaging  methods.  The  use  of  mortar 
elements  is  proven  to  be  optimally  convergent  on  nonmatching  grids,  provided  the  finite  element  space 
used  for  the  interface  variables  consists  of  piecewise  polynomials  of  one  degree  higher  than  the  trace 
along  the  interface  of  the  finite  element  space  used  to  approximate  the  flux  within  the  subdomains  (Ar¬ 
bogast  et  al.,  2000).  Nonetheless,  while  mortar  elements  are  well  known  in  some  application  domains, 
e.g.,  flow  in  porous  media,  they  are  not  widely  employed  for  multiphysics  problems.  Rather,  various 
“geometric”  techniques  are  used  in  most  practical  settings,  especially  in  situations  in  which  one  or  more 
of  the  components  are  solved  with  legacy  “blackbox”  codes.  This  second  approach  is  often  rationalized 
using  a  combination  of  ad  hoc  formal  stability  and/or  accuracy  arguments  combined  with  high  perfor¬ 
mance  computing  expediences.  Moreover,  in  the  situation  in  which  legacy  codes  are  used  to  solve  either 
component,  there  is  little  choice  because  of  the  very  considerable  investment  that  would  be  required  to 
replace  these  codes. 

We  are  not  arguing  for  or  against  either  mortar  elements  or  “geometric”  approaches.  Rather,  we  ad¬ 
dress  two  issues:  (1)  What  effect  do  these  coupling  approaches  have  on  accuracy  of  specified  quantities 
of  interest?  and  (2)  In  each  case,  quantify  the  relative  contributions  of  various  aspects  of  discretization 
to  the  error  in  the  computed  information.  The  tool  we  use  to  address  these  issues  is  an  adjoint-based 
a-posteriori  error  estimate  (Estep  et  al.,  2000;  Becker  &  Rannacher,  2001;  Giles  &  Suli,  2002;  Wheeler 
&  Yotov,  2005;  Estep  et  al.,  2009a;  Hansbro  &  Larson,  2011;  Pencheva  et  al.,  2013).  This  goal-oriented 
estimate  accurately  quantifies  various  contributions  to  the  overall  error.  In  particular,  the  estimate  distin¬ 
guishes  contributions  specifically  arising  from  the  mis-matched  grids  and  the  way  in  which  the  coupled 
information  is  approximated.  We  identify,  through  numerical  examples,  cases  in  which  the  geometric 
projections  are  the  dominant  source  of  error  by  one  to  two  orders  of  magnitude. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  two  introduces  the  differential  equation 
and  the  details  of  the  two  discrete  methods.  Section  three  derives  the  a-posteriori  error  estimate.  Sec¬ 
tion  four  contains  the  numerical  experiments.  Section  five  discusses  computational  logistics  related  to 
iterative  solvers,  and  a  brief  conclusion  is  given  in  section  six. 
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2.  Definition  of  the  problem  and  discretization  methods 

We  define  the  coupled  problem  with  a  common  interface,  then  describe  the  finite  element  and  finite 
volume  discretizations.  We  employ  the  well  known  equivalence  between  finite  volume  methods  and  the 
mixed  finite  element  method  (Russell  &  Wheeler,  1983;  Weiser  &  Wheeler,  1988)  to  recast  everything 
in  the  finite  element  framework.  This  greatly  eases  the  derivation  of  a-posteriori  error  estimates  and 
provides  a  systematic  framework  for  describing  geometric  approaches  to  computing  coupling  values. 

2.1  The  differential  equation 

The  differential  equation  (2.1)-(2.3)  consists  of  a  system  of  second  order  elliptic  partial  differential 
equations  (PDE)  in  two  spatial  dimensions.  The  system  is  posed  on  a  rectangular  domain  T2  consisting 
of  two  nonoverlapping  rectangular  subdomains,  Q/  on  the  left-hand  side  and  Qr  on  the  right-hand  side, 
that  share  a  common  interface  f},  and  whose  union  forms  the  entire  domain,  as  shown  in  Fig.  1.  The 
unit  normal  vector  n  is  defined  to  point  from  left  to  right  on  f/,  and  is  an  outward  pointing  normal 
on  T)  =  \  T/  and  TK  =  r)QK  \  T/.  For  simplicity  of  presentation,  we  assume  Dirichlet  boundary 

conditions  on  dfl,  the  external  boundaries  of  the  domain.  The  results  extend  to  problems  with  Neumann 
conditions  on  part  of  the  boundary  in  a  straightforward  way. 


rL 
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Fig.  1.  Subdomains,  boundaries,  and  definition  of  normal  n  on  the  interface. 


For  a  diffusion  function  a,  split  as  ol  £  W1,°°(f2£)  and  or  £  source  function  /,  split  as 

/ l  G  L2(£2l)  and  fR  £  L2(Qr),  and  boundary  data  g,  split  as  gr  £  H3/2(Il)  and  gR  £  H3I2(Tr),  the 
coupled  system  is 


aL'uL  +  VpL  =  0, 

(x,y)  e  nL: 

<1 

s 

i'¬ 

ll 

> 

(x,y)  £  nL, 

(2.1) 

PL  =  gL , 

(x,y)  £  rLl 

aRiuR  +  VpR  =  0, 

{x,y)  e 

V  -Ur  =/r, 

(x,y)  £  nRl 

(2.2) 

PR  =  gRi 

(x,y)  £TRl 

k  =PL  =  PR , 

(x,y)  £  I], 

(2.3) 

o' 

II 

3 

1 

3 

(■ x,y )  g  17, 
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where  we  assume  that  the  diffusion  matrices,  ol  and  aR,  are  functions  of  space  times  the  identity,  i.e.. 


DL(x,y ) 

0 

DR(x,y) 

0 

0 

DL(x,y)  _ 

0 

DR{x,y)  _ 

(2.4) 


with  Dj  £  W1,00(f2(),  i  —  L,R,  and  min_Z>,(x,y)  ^  Dq  >  0,  so  aj  is  invertible  and  uniformly  coercive 

(x,y)eni 

for  i  =  L,R.  Note  that  we  have  defined  q  as  the  common  interface  pressure  in  (2.3). 

2.2  Mixed  finite  element  mortar  discretization 

The  mortar  finite  element  discretization  was  developed  precisely  for  the  situation  presented  by  dis¬ 
cretization  of  (2.1)-(2.3)  using  two  different  grids  in  the  two  different  subdomains.  We  assume  that 
each  subdomain  is  discretized  by  a  (logically)  rectangular  finite  element  grid.  Lagrange  multipliers  are 
introduced  on  the  interface  boundary  to  provide  a  weak  formulation  of  the  pressure  coupling  conditions. 
Since  the  grids  are  different  on  the  two  sides  of  the  interface,  the  Lagrange  multiplier  space  cannot  be 
the  normal  trace  of  the  velocity  space.  So,  we  introduce  a  mortar  finite  element  space  on  the  inter¬ 
face  (Arbogast  et  al.,  2000;  Bernardi  et  al.,  2005;  Arbogast  et  al.,  2007).  As  shown  in  Arbogast  et  al. 
(2000),  the  method  is  optimally  convergent  and  has  several  other  desirable  convergence  properties  if  the 
boundary  space  has  one  order  higher  approximability  than  the  normal  trace  of  the  velocity  space.  The 
same  order  of  convergence  is  obtained  for  both  continuous  or  discontinuous  piecewise  polynomials  in 
the  mortar  space.  In  our  discretization,  we  choose  the  interface  grid  that  has  one  cell  for  every  two  cells 
in  the  finer  of  the  two  subdomain  grids.  Fig.  2  shows  the  arrangement  for  a  5  x  5  grid  next  to  8  x  8  grid. 
(Note  that  our  convention  is  that  the  finer  grid  is  always  used  in  the  righthand  subdomain.) 


Fig.  2.  Example  grid  shown  separated  into  the  part  on  Qi,  I/,  and  Qr  from  left  to  right. 


We  use  standard  L2  inner  product  notation,  i.e.,  for  functions  F  and  G  defined  on  Q ,  split  as  above, 

{Fi,Gi)=  [  Fi(x,y)  Gj(x,y)dxdy,  i  =  L,R, 

jQj 


and  for  functions  defined  on  the  boundaries,  we  similarly  denote 

(F,G)  r=  [  FGds ,  ;  =  L,I,R. 


The  mixed  finite  element  (mortar)  method  starts  with  the  following  continuous  weak  formulation.  Find 
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Pi  G  Wj  =L2(£2j),  Uj  G  Vi  =  H(div,£2i),  E,  G  A  =  //1//2(f}),  i  =  L,R,  satisfying 

(al'uL,vL)-{pL,V  -vL)  +  (^n-vL)ri  =  ~(gL,n  ■  vL)rL, 

(V-uLlwL)  =  {fL,wL), 

( aR'uR,vR )  -  (pR,  V  •  vR)  -  (< %,n  vR)r,  =  -{gR,n-  vR)rR, 

(V-uR,wR)  =  (fR,wR),  (2.5) 

{n-{uL-uR),v)r,  =0, 
for  all  ( Wj,Vj,V )  G  ( Wi:Vi,A ),  i  =  L, R. 

To  discretize,  we  use  the  lowest  order  Raviart-Thomas  finite  element  space  (RT 0),  in  which  the 
discrete  scalar  unknown  ph  is  approximated  as  a  constant  over  each  cell,  and  the  components  of  the 
discrete  vector  ub  are  approximated  by  functions  that  are  piecewise  linear  in  one  spatial  dimension  and 
constant  in  the  other  (Bernardi  et  al. ,  2005;  Estep  et  al. ,  2009a).  The  discrete  interface  unknown,  qh ,  is 
represented  by  piecewise  discontinuous  linears  on  the  interface  grid  cells  (Arbogast  et  al.,  2000,  2007). 
The  test  functions  in  the  discretization  of  the  weak  formulation  of  (2.5)  corresponding  to  w,  v,  and  v 
are  restricted  to  these  same  spaces.  To  be  precise,  for  a  finite  element  partition  A  of  [a,/?],  and  for 
r  =  0, 1,2, ....  q  =  —1,0, 1, ...,  we  define  the  piecewise  polynomial  space 

^{A)  =  {vG  Cq([a,b\) : 

v  is  a  polynomial  of  degree  C  r  on  each  subinterval  of  A  } . 

When  q  =  —  1  the  functions  are  discontinuous.  The  space  of  continuous  piecewise  bilinear  functions  is 
the  tensor  product  ( Ax )  (g>  (Ay).  The  RT0  discrete  spaces  are 

W?  =.//°l{AX,i)®^°i(.Ay4),  i  =  L.R. 

V)  =  (AXj)  C>  1  (Ajy)]  x  [,-#° j  (AXti)  <g>  (AVj,-)] ,  i  =L,R, 

A1'  =  [  ( Aj-j ) . 

The  mixed  finite  element  (mortar)  method  reads:  Compute  pb  G  Wj\  ub  G  V1-,  c,h  G  Ah,  i  =  L,R, 
satisfying 

(aL 1  “l »  Vl)  -  {phL,  V  •  vL)  +  (? h ,  n  ■  vL)r,  =  - (gL , «  •  vL)rL , 

=  (fL,wL), 

(aR 1  uhR,  vR)  -  (pR,  V  ■  vR)  -  (|  ■ * ,  n  ■  vff )r,  =  -(gR,n  ■  vR)r,t , 

(V-ur,wr)  =  (fR,wR),  (2.6) 

(n-(u'[-uhR),v)ri=  0, 

for  all  (w;,  V;,  v)  G  (Wb ,Vb ,Ah),  i  =  L.R.  This  yields  a  discrete  system  of  the  form 


1 

bs 

o 

o 

’  UL  ' 

'  -Dl  ' 

Btl  0  0  0  0 

Pl 

Fl 

0  0  Mr  -Br  Cr 

uhR 

= 

-Dr 

0  0  Btr  0  0 

p'r 

Fr 

o 

o 

o 

- 1 

[  zh  \ 

0 

where  we  abuse  notation  to  let  ub,  pb,  and  qh  denote  the  vector  of  nodal  values  for  the  finite  element 
functions. 
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2.3  Geometrically  coupled  finite  volume  discretization 

The  standard  formulation  of  the  finite  volume  method  eschews  a  variational  formulation  of  the  problem, 
so  there  is  no  natural  description  of  a  weak  imposition  of  the  coupling  conditions  in  that  formulation. 
Moreover,  the  standard  finite  volume  method  provides  approximation  values  of  p  only  at  cell  cen¬ 
ters  while  approximate  values  for  u  along  cell  boundaries  are  obtained  by  differencing  the  p  values. 
These  characteristics  motivate  the  use  of  “geometric”  coupling  techniques  that  employ  a  combination 
of  extrapolation  and  averaging  to  provide  coupling  values  of  both  unknowns  along  the  interface.  The 
motivation  for  this  approach  is  reinforced  in  the  context  of  iterative  solution  of  the  coupled  problems, 
where  well  posed  problems  are  created  on  each  subdomain  using  interface  boundary  conditions  ob¬ 
tained  from  the  other  subdomain.  In  this  approach,  it  is  necessary  to  couple  the  coarser  side  using  state 
values  extrapolated  from  the  finer  side  solution,  while  the  finer  side  must  be  coupled  to  flux  values, 
which  are  themselves  differences  of  state  values,  extrapolated  from  the  coarser  solution.  Reversing  this 
arrangement  can  lead  to  a  singular  system. 


Fig.  3.  Extrapolation  to  the  interface.  Left:  Neumann 
values  on  the  interface  are  computed  by  linear  extrapola¬ 
tion  of  the  last  two  available  flux  values,  which  are  dif¬ 
ferences  of  state  values.  Right:  Dirichlet  values  on  the 
interface  are  computed  by  linear  extrapolation  of  the  last 
two  available  state  values. 


Fig.  4.  Averaging  or  broadcasting  of  extrapolated  values.  Left:  In 
the  case  of  constant  extrapolation,  the  last  available  state  or  flux 
value  is  simply  used  as  the  interface  value.  Right:  Weighted  av¬ 
eraging  of  state  and  flux  values  when  cell  widths  do  not  share  an 
integer  ratio. 


To  obtain  values  on  the  interface,  we  employ  either  linear  or  constant  extrapolation.  We  illustrate 
linear  extrapolation  in  Fig.  3.  We  compute  the  extrapolated  values  by  computing  a  linear  or  constant 
interpolant,  which  is  then  evaluated  at  the  interface  boundary.  We  denote  the  extrapolated  values  using 
the  operators  Pr~>l(Pr)  rrncl  PL^R(pbL).  When  the  cells  on  either  side  of  the  interface  do  not  match,  then 
weighted  averaging  and  “broadcasting”  schemes  are  used  to  generate  values.  In  Fig.  4,  we  illustrate  the 
averaging  and  broadcasting  schemes  when  two  cells  on  the  right  match  one  cell  on  the  left.  The  state 
values  at  the  two  circle  locations  are  averaged  and  used  at  the  square  location.  The  flux  value  at  the 
square  location  is  “broadcast”  to  both  of  the  circle  locations.  When  the  cell  widths  on  the  coarse  and 
fine  side  of  the  interface  do  not  share  an  integer  ratio,  then  a  suitable  averaging  of  values  is  used.  For 
example,  in  the  2  cells  next  to  3  cells  arrangement  pictured  in  Fig.  4,  the  state  value  at  location  D  is  set 
equal  to  |  the  state  value  at  location  A  plus  i  the  state  value  at  location  B.  The  flux  value  at  location  A 
is  set  equal  to  the  flux  value  at  location  D,  while  the  flux  value  used  at  location  B  is  set  equal  to  half  the 
flux  value  at  D  plus  half  the  flux  value  at  E. 

We  formulate  the  finite  volume  method  as  an  RTO  mixed  finite  element  method  employing  a  special 
quadrature  formula,  following  Russell  &  Wheeler  (1983);  Weiser  &  Wheeler  (1988).  This  provides 
a  foundation  for  deriving  an  a-posteriori  error  analysis  for  the  finite  volume  scheme,  see  Estep  et  al. 
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(2009a).  The  version  of  (2.6)  equivalent  to  a  finite  volume  method  reads:  Compute  p1'  G  W-1,  i/f  €  Vf, 
tqh  €  Ah,  i  =  L,R,  satisfying 

(' alluhL,vL)Mj  ~  vL)  +  (PR^L{pR),n  •  vL)r,  =  ~{gL,n  •  vL)rL.M, 

(V  •  UhL,wL)  =  (fL,wL), 

(aRlUR,VR)M, t  -  (Pr,V-vr)  vR)r,  =  -( gR,n  ■  vR)rR,M, 

(V  •  Ur,wr)  =  (/r,wr),  (2.8) 

({Pl^r(p'L)  -  n  ■  4),v)r,  =  0, 

for  all  (w;.  V;.  V)  G  ( Wj1  ,V'j  ,Ah),  i  =  L,R .  Here  we  employ  the  approximate  inner  product 

(u'\v)m,t  =  {Ux,Vx)Tx,My  +  (Uy,Vy)Mx,Ty, 


where  ;V/: . :  and  77  denote  the  midpoint  and  trapezoidal  quadrature  rules  in  the  x  and  y  directions  as 
indicated,  while  (• ,-)rhM  denotes  the  midpoint  rule  for  i  =  L,  R.  Note  the  quadrature  formulas  are  applied 
internally  on  each  cell,  so  potential  discontinuites  in  a  and  f  across  f}  cause  no  difficulty. 

This  yields  a  discrete  system  of  the  form 


'  Ml  -Bl  0  Qd  0 

’  U'L  ' 

'  -Dl  ' 

Btl  0  0  0  0 

Pl 

Fl 

0  0  Mr  -Br  Cr 

UR 

= 

-Dr 

0  0  Btr  0  0 

Pr 

Fr 

o 

to 

a 

& 

o 

o 

l  \ 

0 

which  should  be  compared  to  (2.7). 

It  is  possible  to  eliminate  the  unknowns  u1-,  i  =  L,R ,  and  q1' ,  to  reduce  (2.9)  to  a  system  for  p1’  of  the 
form 


'  AL 

Cd  ' 

'  FL  ' 

Cn 

Ar 

[p'r  \ 

.  pR  - 

(2.10) 


The  averaging  and  broadcasting  are  incorporated  into  the  “coupling  Dirichlet”  and  “coupling  Neumann” 
matrices  Cd  and  C,y-  This  is  the  same  system  that  is  constructed  by  using  a  finite  volume  approach 
directly. 

We  have  verified  through  numerical  experiments  that  the  p  component  of  the  solution  of  (2.9)  is 
identical  to  the  solution  of  (2.10).  Furthermore,  the  u  component  of  the  solution  of  (2.9)  is  identical  to 
the  u  values  obtained  by  differencing  the  solution  of  (2.10)  to  approximate  V p  at  the  cell  boundaries 
and  evaluating  the  diffusivity  at  the  cell  boundaries.  The  q  component  of  the  solution  of  (2.9)  has  no 
counterpart  in  the  solution  of  (2.10). 


3.  A-posteriori  error  analysis 

Our  goal  is  to  derive  an  a-posteriori  error  estimate  for  the  quantity  of  interest 

(euL-,¥uL)  +  (ePL>  Vpl)  +  (e“R'¥uR)  +  (ePR’¥pR)  +  (e% )  Vz)r„  (3.1) 

where  l/ru, .  \j/pL .  \jfUR .  i/r/)R ,  and  are  given  L2  functions  and  e(.;  denotes  the  errors  in  the  corresponding 
variables.  We  define  the  generalized  Green’s  function  corresponding  to  these  functionals  using  the 
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adjoint  problem 


a 

t-1 

-©- 

t-i 

1 

<1 

II 

R 

C~” 

on  12L, 

-V- <I>l  =  Vpl 

on  12L, 

(3.2) 

o 

II 

Co 

on  11, 

=  v^U/( 

on 

-V-<t>R  =  Vpr 

on  12^, 

(3.3) 

o 

II 

on  11, 

(3  =  Cl  =  Cr 

on  1), 

(3.4) 

"■(Ql-Qr)  =  Vs 

on  1/. 

The  a-posteriori  error  estimates  explicitly  depend  on  <j>L,( £l,<I>r,  and  £r. 


3. 1  Estimate  for  mortar  mixed  finite  element  method 

We  first  derive  an  estimate  for  the  mortar  finite  element  method  assuming  all  integrals  in  the  weak 
formulation  are  computed  exactly.  We  begin  by  substituting  (3.2)— (3.4)  for  the  various  i//’s  in  (3.1)  and 
applying  the  divergence  theorem, 

(euLiWuL)  +  (ePL>VpL)  +  V%)  +  (ePRiWpR)  +  (eS  >  Vs)n  (3.5) 

=  ( eUL,all  <t>L)  +  (V-eUL,CL)-{n- eUL1fi)r,  -  (ePL,  V  •  <j>L) 

(euR^R1  <I>r)  +  (y  ■  eUR,Cii)  +  {n  ■  eUR, P)r,  -  (ePRy  ■  <I>R) 

+(e$,n- (Ql-Qr))^- 

Expanding  on  the  right  and  subtracting 

(aL  V,  <t>L)  -  (pL,V-  4>l)  +  (Z »« '  <h)n  +  (gL,n-  $L)rL 

+(aR1uRl<l>R)  -  (pR,V-  $R)-(%,n-  <t>R}ri  +  (, gR,n  ■  <t>R)rR 

+  (V-Mr,Cr)-(/RiCr) 

-(n-(uL-uR),P)r,  =0, 

obtained  by  substituting  the  adjoint  solution  as  test  functions  into  the  forward  weak  form  (2.5),  gives 

(e«L  i  VuL  )  +  (ePL  >  Vpl  )  +  (eUR  >  )  +  (ePR  >  VpR  )  +  (eS )  )f? 

=  ~(aZ' 1«i  0l)  +  (p'L  v  •  0l)  -  •  § l)Il  -  ■  0 i)r; 

+(/L,CL)-(V-«i,CL) 

-(aR'uhR,$R)  +  (4, V ■  -  <g*,»  ■  tf>/(.)rR  +  <€*,»•  0«)r; 

+(/*>&)-(?•«*,&) 

+(n-(uhL-uhR),p)rr 


(3.6) 
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We  rewrite  this  as 

(euL  i  VuL )  +  (ePL  >  Ypl  )  +  (e»R  •  VUR )  +  iePR  >  YpR  )  +  )r;  (3.7) 

=  (^«z.  >  0l)  +  (-^pl  >  Cl)  +  (•*««  ■  Qr)  +  (rpr  >  Cr)  +  (R^  P  )r, , 

wherein  the  residuals  are  given  by 

Ri<L  =  ~aLluL-vPL>  RuR  =  -aR1uhR-VPR > 

RPL=fL-V  -uhL,  RPR=fR-VuhR ,  R^=n-{uhL-uhR). 

Note  that  the  divergence  theorem  implies 

(rul,<I>l)  =  -(aLlu'[, $L)  +  {phL,  V •  0t)  -  {phL,n ■  <j>L)dQL 

=  -i.alluhLJ<\>L)  +  (^,  V  ■  0L)  -  (gL,n  •  0L)rL  -  (<^,n  •  0L)r;- 

Also  note  that  [5  =  =  £r  for  the  continuous  adjoint  solution,  but  /j  is  distinct  from  £/,  and  Ca1  for  the 

discrete  solution. 

Next,  we  use  Galerkin  orthogonality.  We  introduce  projection  operators  that  map  into  the  finite 
element  space  of  the  discrete  forward  solution: 

Pt-.L\nL)^w'l,  Pr  :  L2(Qr)  — »  Wr, 

n'l :  L2(Ql)  -a  Vi,  IIr  :  L2(£>*)  ->  V* ,  Z/!  :  L2(r7)  -A  Aft. 

The  actual  choice  of  projection  is  immaterial  for  the  estimate.  In  practice,  we  employ  a  combination  of 
restriction  and  averaging.  Without  quadrature,  Galerkin  orthogonality  for  (2.6)  is  expressed  as 

(RuL,nfoL)  +  (rpl,PlCl)  +  (RUR,n^<t>R)  +  (rpr,p^r)  +  (R4,zhp)n  =  o, 

and  subtracting  gives  the  following  result. 

Theorem  3.1  The  errors  for  the  mixed  finite  element  method  (2.6)  without  quadrature  satisfy 

(ePL>  Vpl)  +  (euL'VuL)  +  (epRiYpR)  +  (euRiVuR)  +  (3.8) 

=  (rUl  ,+L-nfrL)  +  (RPL,a-p££) 

+  (RUR,<l>R-n^R)  +  (rpr  ,  c R  -  pH  &)  +  <*? ,  j3  -  zh p )r, , 

wherein  the  quantities  on  the  right-hand  side  are  computable  provided  the  true  adjoint  solution  is  avail¬ 
able. 

In  practice,  we  employ  a  numerical  solution  of  the  adjoint  problem.  To  emphasize  this,  we  state  the 
following  corollary  that  involves  numerical  adjoint  quantities. 

Corollary  3.1  Provided  that  the  projection  operators  p£,  P% ,  IT*,  17^,  and  Z/;  are  bounded  in  L2,  the 
errors  for  the  mixed  finite  element  method  (2.6)  without  quadrature  can  be  estimated  as 

(ePii  Ypl)  +  (euL,VuL)  +  (epR’YpR)  +  (euR,YUR)  +  (e£, >  W%)r,  (3.9) 

« i*uL, 4>t  nWl)  +  (RPL, g  - PHCH) 

+  (rur  ,*hR-nfrhR)  +  ( rpr  ,  C H-pH  CH  )  +  <*  5 ,  Ph -  zi,ph)r, , 

for  numerical  solutions  <j)L  ss  </)(',  sa  £  * ,  <f>R  ss  0^,  Cr  <^',  and  P  ~  Ph-  In  this  approximation,  the 
errors  are  to  be  measured  in  the  L2-norm. 
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The  proof  follows  from  the  triangle  inequality  and  the  definition  of  the  operator  norm.  That  is,  the 
absolute  value  of  the  difference  between  the  right-hand  sides  of  (3.8)  and  (3.9)  is  bounded  by 

(i+iin*ii)ii«Bj|2||^-#*i|2+(i+ii^ii)ii/?pj|2|ia-tfii2 
+  (i  +  II^IDII^Ibll^  -  ^lb  +  (i  +  K!II)PpJ2|ICk  -  A 

+  (1  +  ||Z/l||)|!^||2,r/||/3-j3"!|2,r;. 

In  order  to  obtain  accurate  estimates,  the  numerical  adjoint  solutions  must  be  sufficiently  accurate. 
Generally  this  is  satisfied  by  solving  the  adjoint  problems  either  using  a  higher  order  numerical  method 
or  using  a  mesh  sufficiently  refined  from  the  one  used  for  the  forward  discretization.  In  the  context  of 
finite  volume  discretizations,  the  second  approach  is  generally  easier  to  implement.  In  our  numerical 
examples  we  use  a  finer  grid,  and  the  accuracy  of  this  approach  is  illustrated  in  section  4.1. 

3.2  Estimate  for  finite  volume  methods  using  geometric  coupling 

3.2.1  The  effect  of  quadrature.  We  first  derive  an  estimate  for  the  mixed  finite  element  method  (2.6) 
with  quadrature,  which  can  be  applied,  say,  if  f,  gj,  and  a,  are  continuous  in  Q,.  i  =  L.  R.  With  quadra¬ 
ture,  Galerkin  orthogonality  is  expressed  as 

(Ru,M[$l)q  +  (KPl-PI%l)q 

+  {RurMk$r)q  +  (RPr,PrC,r)q  +  (R^,Zh[T}(QTl  =  0, 

where  we  use  the  subscript  Q  to  denote  the  approximate  inner  product  using  quadrature.  It  is  impor¬ 
tant  to  distinguish  residuals  associated  with  approximating  the  solution  spaces  using  finite  dimensional 
polynomial  spaces  from  residuals  associated  with  approximating  the  integrals  defining  the  variational 
formulation.  We  rewrite  Galerkin  orthogonality  as 

{RUL,n£*L)  +  (Rpl^Cl)  +  (*UR,nfrR)  +  (rpr,A)  +  (R^,zh .fi)n 

-  QEUL(n£+L)  -  qepl(p^l)  -  QEuR(nfrR )  -  Qepr(p^r)  -  qe :s  (zhp)  =  o, 

with 

QEUL(n^L)  =  (RULMj‘<t>L)  -  {RUL,ritfL)Q, 

QEpl(P^l)  =  (Rpl,p£Cl)  -  (Rpl,p£Zl)q, 

QEUR(nfoR)  =  (RUR,n'M  -  (RUR,n^R)Q, 

QePr(PrCr)  =  (RPr,PrCr)  -  (Rpr,Pr£r)q, 

QE.(Zhp)  =  (R^Zh P)n  -  (R^Zh p)QJ}. 

This  gives  the  following  a-posteriori  estimate  for  the  mixed  finite  element  method  with  quadrature. 

THEOREM  3.2  If  fi,  gi ,  and  a,  are  continuous  in  Q,.  i  =  L.  R,  then  the  errors  for  the  mixed  finite  element 
method  (2.6)  with  quadrature  satisfy 

( ePL  >  Vpl)  +  (e«L  >  Vul  )  +  ( epR  >  Vpr  )  +  (e«R  >Vas)+  (e$  >  Wf,  )n 
=  (*uL ,  0L  -  nfrL)  +  ( rPl  ,  Cl  -  p£Zl) 

+  (rUr  ,  <I>R  -  n^R )  +  (RPR ,  c«  -  p£  Cr)  +  (R^P-  zhp)r, 

+  QEUL(n^L)  +  QEPl(Pp  Cl)  +  QEUR  (n^R)  +  QEPr  (Pk  C«)  +  QEt;  (Z/!j3 ) . 


(3.10) 
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Note  that  in  the  case  of  using  the  RTO  finite  element  space  and  the  midpoint-trapezoidal  quadrature 
rules  discussed  above,  the  mixed  finite  element  method  reduces  to  the  finite  volume  method  (Russell  & 
Wheeler,  1983;  Weiser  &  Wheeler,  1988;  Estep  et  al.,  2009a),  and  some  of  the  quadrature  error  terms 
are  zero.  These  terms  are  included  for  generality,  so  that  (3.10)  is  valid  for  other  combinations  of  finite 
element  spaces  and  quadratures. 

Note  that  in  practice,  we  implement  the  obvious  analog  of  Corollary  3.1,  which  now  requires  suffi¬ 
cient  smoothness  of  the  solution  to  obtain  sufficiently  accurate  quadrature  approximations. 

3.2.2  The  effect  of  geometric  coupling.  For  the  geometric  coupling  (2.8),  the  Galerkin  orthogonality 
becomes 

(RUL,nfrL)Q  -  (pR^L(phR)-^h,n-n^L)Qri  +  {rpl,p£Cl)q 

+  (RuR,njfrR)Q  +  0 Rpr,p'Kr)q  +  (RsMan  -  (» •  uhL-pL^R(PhL),zh j3)e,r;  =  o. 


Defining 


2*UL{nfrL)  =  {RUL,n}ML)-{RUL,nltL)Q 

+  (PR^L(l 4)  -  Sh,nfrL)Qri  -  (pR^L(phR)  -  •  n^L)r„ 

i)  =  (R^ZhP'fr,  -  (R^,ZhP)Qjj 

+  (nL-uhL-  PL^R {phL),Zhfi)Qf,-  (n-uhL-  PL^R (phL), Zh  fi)rn 

and  arguing  as  above  gives  the  following  result. 

Theorem  3.3  If  fj,  gi,  and  a,  are  continuous  in  Q,,  i  =  L,R ,  then  the  error  for  the  mixed  geometric 
finite  volume  method  (2.8)  satisfies 

(ePi>  Vpl)  +  (e«L’  VV)  +  ( ePR’V/PR )  +  ieuR,VUR)  +  (e^Y%)rI  (3.11) 

=  (R*l,*l  -  n^L)  +  (Pr^l(p'r)  -^\n-  U^L)ri  +  (Rpl,Cl  - /£&) 

+  (rur  ,<i>R-n^R)  +  (RPR ,  c«  -  P'1  C R) 

+  (R$, P  - zhP ) r,  +  (« •  - Pl^r ( Pl ) ,ZhP) r, 

+  ^UL(n^L)  +  qePl{p£  Cl)  +  QeUr  (nfoR)  +  qePr  (/*  &  )  +  ^  (zh  p ) . 

Note  that  in  practice,  we  implement  the  obvious  analog  of  Corollary  3.1,  assuming  again  sufficient 
smoothness  of  the  solution  to  obtain  sufficiently  accurate  quadrature  approximations. 

4.  Numerical  investigations 

In  this  section,  we  use  the  a-posteriori  error  estimates  to  investigate  in  detail  the  accuracy  of  the  two 
approaches  to  coupling.  For  all  of  the  investigations,  the  coarser  subdomain  Eli  is  given  by  x  £  [—1,1] 
and  y  £  [—2,0],  the  finer  subdomain  QR  is  given  by  x  £  [—1,1]  and  y  £  [0,2]  (see  Fig.  1),  and  the 
interface  f)  is  located  along  y  =  0.  (Note  that  here  the  bottom  subdomain  is  considered  as  being  “left” 
and  the  top  one  is  “right,”  in  conformance  to  our  convention  as  to  the  finer  subdomain.)  The  grids  are 
reported  as  nR  x  mR  for  the  left  domain  and  nR  x  mR  for  the  right  domain,  where  uia  corresponds  to 
the  number  of  cells  in  the  x-direction  (which  is  also  the  number  of  cells  along  the  interface),  and  m,.  . 
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corresponds  to  the  number  of  cells  in  the  ^-direction.  The  boundary  conditions  for  all  tests  are  Dirichlet. 
To  avoid  issues  arising  from  iterative  solution  of  the  discrete  system,  we  employ  direct  methods  to  find 
the  approximate  solution  to  within  machine  precision. 

The  quantity  of  interest  being  sought  is  specified  by  giving  the  adjoint  problem  data  \j/Ux,  i \rUy,  y/p, 
and  y/= .  The  adjoint  problem  is  solved  using  the  same  RTO  mixed  finite  element  method,  but  on  a  grid 
that  is  significantly  finer  than  that  of  the  forward  problem,  so  that  the  discretization  error  associated  with 
the  adjoint  solution  has  no  significant  effect  on  the  results. 

The  functions  chosen  for  the  source,  diffusivity,  and  adjoint  data  are  either  constants  or  Gaussian 
functions  of  the  form 


e-(y~br 


+  K, 


which  gives  a  localized  “ridge”  centered  at  y  =  b.  In  the  case  of  the  adjoint  data,  the  Gaussian  or  constant 
function  being  used  is  normalized  so  that  the  area  under  if/  is  equal  to  one.  The  parameter  K  is  non  zero 
only  in  the  case  of  diffusivity,  where  this  constant  is  added  to  the  Gaussian  to  prevent  the  diffusivity 
from  approaching  zero  anywhere  in  the  domain. 

In  the  tests,  we  report  values  for  the  terms  in  (3.10)  and  (3.11)  that  are  non  zero.  For  both  the  mixed 
finite  element  and  geometric  finite  volume  methods  the  following  five  terms  are  included: 


MFEi 

or 

GFV\  =  (Ru,  ,<l>b  -  Flj'Q1}) , 

mfe2 

or 

GFv2  =  (RUR,<t>hR-n *thR), 

MFEi 

or 

GFV3  =  (RPL,tf-Pbtt), 

MFEs, 

or 

GFV4  =  (RPR,tf-Pktf), 

mfe5 

or 

GFV5  =  {R^ph-Zhph)rr 

In  the  geometric  finite  volume  case,  we  add  two  additional  terms  relating  to  the  geometric  projections 
and  two  additional  quadrature  terms: 


GFV6  =  {Pr^l(A)  -  ? \n  ■  n^'Dr,, 
GFVi  =  ( n  ■  uhl-Pl^R(phL),Zhfih)rn 
GFVz=£gUL{nhLthL), 

GFV9  =  QEUR(T1 hRthR). 


We  note  that  the  first  five  expressions,  common  to  both  MFE  and  GFV,  are  often  similar  in  size.  As  a 
gross  measure  of  the  effect  of  geometric  projection  and  of  the  use  of  quadrature,  we  also  report  the  two 
ratios 


Lie  I  GFV,  | 


Lit  |  GFVi  | 


ratiOqUad 


Lis  I  GFVj  | 
Lli  I  GFV,  I ' 


We  present  three  examples  chosen  to  show  the  spectrum  of  possibilities  in  terms  of  the  perfor¬ 
mance  of  the  methods.  Test  case  1  is  an  “easy”  case  in  which  the  geometric  method  performs  relatively 
well.  Test  case  2  has  a  narrow  and  severe  dip  in  diffusivity  located  along  the  subdomain  boundary  and 
consquently  the  geometric  method  performs  poorly,  as  might  be  expected  for  a  problem  in  which  the 
solution  changes  rapidly  near  the  interface.  Test  case  3  is  based  loosely  on  a  real  world  fusion  problem 
and  demonstrates  one  of  our  main  conclusions,  which  is  that  the  geometric  method  can  perform  poorly 
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even  when  the  solution  is  smooth  near  the  interface.  Note  that  the  behavior  of  the  diffusivity  function 
approaching  the  common  interface  on  both  sides  has  more  impact  on  the  accuracy  of  coupling  than  a 
discontinuity  in  the  diffusivity  across  the  interface. 

4.1  Verification  of  a-posteriori  estimate  accuracy 

We  begin  with  a  problem  for  which  we  have  manufactured  the  known  solution 

P{x,y)  =c°s(y)cos(^p).  (4.1) 

The  diffusivity  a  is  equal  to  one  everywhere.  The  other  solution  components,  the  source  term  /,  and  the 
boundary  values  g  for  the  problem  follow  from  (4.1).  Since  we  know  the  true  solution,  we  can  compute 
the  exact  error  terms  (e,  yr)  on  the  left  in  (3.10)  and  (3.11)  directly  and  then  compare  to  estimates  of 
the  quantities  on  the  right  computed  using  a  numerical  solution  to  the  adjoint  problem.  In  this  situation, 
the  most  important  issue  for  the  accuracy  of  the  estimates  is  the  accuracy  of  the  approximate  adjoint 
solutions.  As  the  grid  for  the  adjoint  problem  is  refined,  the  estimates  become  more  accurate.  That  is, 
using  the  approximation  to  the  adjoint  problem,  the  estimated  quantities  Y.MFEi  or  }_ CFV,  becomes 
closer  to  their  true  value,  the  error  in  the  quantity  of  interest  MFE  £(e,  yr)  or  GFV  }_(e.  y/).  Tables  1 
and  2  show  this  using  coarse  and  fine  forward  solutions. 


Table  1.  The  forward  problem  with  solution  (4.1)  is  run  at  10  x  10  next  to  16  x  16.  The  adjoint  problem  is  run  at  several  grids 
to  show  how  the  sum  of  terms  approaches  the  direct  calculation  of  (e,  if/).  The  adjoint  data  components  i \ru  and  i f/p  are  constant 
everywhere  and  i//^  =  0. 


adj.  grid 

MFEX>,  y/) 

ZMFE, 

ratio 

GFV  £(e,  y) 

£  GFV i 

ratio 

20x20  :  32x32 

1.96£— 3 

1.47£-3 

.749 

—  1.00E—3 

— 1.50£-3 

1.49 

40x40  :  64x64 

1.96E-3 

1.84£— 3 

.937 

—  1.00E—3 

—  1.13E-3 

1.12 

80x80  :  128x128 

1.96E-3 

1.93E-3 

.984 

—  1.00E—3 

—  1.03E— 3 

1.03 

160x160  :  256x256 

1.96E-3 

1.96£— 3 

.996 

—  1.00E—3 

— 1.01£— 3 

1.01 

Table  2.  The  forward  problem  with  solution  (4.1)  is  run  at  40  x  40  next  to  64  x  64.  The  adjoint  problem  is  run  at  several  grids 
to  show  how  the  sum  of  terms  approaches  the  direct  calculation  of  ( e ,  i f/).  The  adjoint  data  components  \f/u  and  t ffp  are  constant 
everywhere  and  =  0. 


adj.  grid 

MFE£(e,  Y) 

ZMFEi 

ratio 

GFV  £(e,  Y) 

Y.GFV, 

ratio 

80x80  :  128x128 

1.23£-4 

9.23E-5 

.750 

—1.00E—5 

— 1.01£— 4 

1.44 

160x160:  256x256 

1. 23£-4 

1.15£-4 

.937 

—1.00E—5 

—1.11E—5 

1.11 

4.2  Convergence 

To  compare  the  accuracy  of  the  various  approximations,  we  use  the  2-norms 


Kll2  =  ^K-^)2,  IM2  =  /£($-Sa)2- 
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We  use  the  manufactured  solution  from  the  previous  section  (a  =  1  and  p  is  given  by  (4.1)).  We  compare 
the  2-norm  errors  of  the  finite  element  and  geometric  finite  volume  methods  on  a  sequence  of  grids  in 
order  to  asses  the  convergence  rate.  The  coarsest  grid  is  10  x  10  next  to  16  x  16,  and  the  number  of  cells 
in  each  dimension  is  doubled  with  each  refinement. 

The  results  in  Tables  3-6  show  that  the  convergence  rate  for  the  geometric  finite  volume  deteriorates 
for  the  ux,  uy,  and  |  components  when  the  number  of  cells  along  the  fine  side  of  the  interface  is  not  an 
integer  multiple  of  the  number  of  cells  along  the  coarse  side  of  the  interface.  When  the  test  is  repeated 
with  a  grid  starting  at  8  x  8  next  to  16  x  16,  the  convergence  rates  for  the  two  methods  are  equal.  The 
first  order  convergence  of  p  and  u  for  the  MFE  is  to  be  expected  (Arbogast  et  ai,  2000). 


Table  3.  Convergence  of  solution  component  p,  indicating  a  rate  of  0(h). 


grid 

MFE  ||<?p|| 

MFE  ratio 

GFV  ||e„|| 

GFV  ratio 

10x10:  16x16 

1.20E-01 

N/A 

1.20E-01 

N/A 

20x20 :  32x32 

5.98E-02 

2.00 

5.98E-02 

2.00 

40x40  :  64x64 

2.99E-02 

2.00 

2.99E-02 

2.00 

80x80  :  128x128 

1.49E-02 

2.00 

1.49E-02 

2.00 

160x160 :  256x256 

7.47E-03 

2.00 

7.47E-03 

2.00 

Table  4.  Convergence  of  solution  component  ux,  indicating  a  rate  of  about  0(h). 


grid 

MFE  IKH 

MFE  ratio 

GFV  IK  || 

GFV  ratio 

10x10  :  16x16 

8.49E-02 

N/A 

8.63E-02 

N/A 

20x20  :  32x32 

4.21E-02 

2.02 

4.26E-02 

2.02 

40x40  :  64x64 

2.10E-02 

2.00 

2.14E-02 

1.99 

80x80 :  128x128 

1.05E-02 

2.00 

1.09E-02 

1.97 

160x160  :  256x256 

5.25E-03 

2.00 

5.63E-03 

1.93 

Table  5.  Convergence  of  solution  component  uy.  indicating  a  rate  of  0(h)  for  MFE  but  less  for  GFV. 


grid 

MFE  IK  || 

MFE  ratio 

GFV  K II 

GFV  ratio 

10x10  :  16x16 

8.41E-02 

N/A 

8.59E-02 

N/A 

20x20  :  32x32 

4.20E-02 

2.00 

4.39E-02 

1.96 

40x40  :  64x64 

2.10E-02 

2.00 

2.29E-02 

1.92 

80x80 :  128x128 

1.05E-02 

2.00 

1.23E-02 

1.86 

160x160  :  256x256 

5.25E-03 

2.00 

6.94E-03 

1.77 

Table  6.  Convergence  of  solution  component  £,  indicating  a  rate  of  0(h~)  for  MFE  but  only  O(h)  for  GFV. 


grid 

MFE  Kl 

MFE  ratio 

GFV  ||e,  || 

GFV  ratio 

10x10:  16x16 

7.53e-03 

N/A 

6.11e-03 

N/A 

20x20 :  32x32 

1.89e-03 

3.99 

1.79e-03 

3.42 

40x40  :  64x64 

4.72e-04 

4.00 

6.37e-04 

2.80 

80x80  :  128x128 

1.18e-04 

4.00 

2.77e-04 

2.30 

160x160:  256x256 

2.95e-05 

4.00 

1.33e-04 

2.09 
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4.3  Test  Case  1 

In  the  next  problem,  we  explore  accuracy  for  a  solution  that  is  not  changing  rapidly  near  the  interface. 
We  find  that  the  use  of  geometric  projections  does  not  lead  to  significant  effects  on  accuracy.  We  let  the 
diffusivity  a  be  one  in  both  £2/  and  £2 r  and  use  the  manufactured  solution  given  by  (4.1).  The  grid  for 
the  forward  problem  is  20  x  20  next  to  32  x  32.  The  adjoint  grid  is  80  x  80  next  to  128  x  128,  and  the 
adjoint  data  is  a  nonzero  constant  for  y/Ux,  i jfUy,  and  y/p,  while  y/^  =  0. 

We  list  the  error  contributions  in  Table  7.  For  the  geometric  approach,  we  list  results  for  both 
constant  and  linear  extrapolation.  The  results  show  that  the  projection  error  for  linear  extrapolation  is 
only  about  one  quarter  of  the  residual  error,  while  the  projection  error  for  constant  extrapolation  is  much 
larger.  Fig.  5  shows  the  solution  components  for  the  finite  element  case.  The  geometric  finite  volume 
solutions  are  very  similar.  Fig.  6  shows  the  adjoint  solution  components. 


Table  7.  Error  terms  for  Case  1.  The  forward  grid  is  20  x  20  next  to  32  x  32.  The  adjoint  grid  is  80  x  80  next  to  128  x  128. 


term 

MFE 

GFV  (linear) 

GFV  ( constant ) 

i 

(RUL,rL-n£K) 

—  1 . 6i? — 4 

—  1.675—4 

—  1.575—4 

2 

(RUli.K-n£K) 

—6.1E—5 

—6.175—5 

—6.175—5 

3 

(Rputt-nzi) 

4.975—4 

4.975—4 

4.975—4 

4 

1 .975—4 

1.975—4 

1.975—4 

5 

(Rz,Ph-Zi,P")rI 

4.275—8 

—  1.375—6 

1.075—5 

6 

(PR-*Ltfh)-Zk,n-nkL*hL)r, 

N/A 

2.075—4 

1.775—4 

7 

(«■«'/ 

N/A 

2.275—5 

3.975—3 

8 

N/A 

—7.175—4 

—7.075—4 

9 

QEuAn'M) 

N/A 

—2.875—4 

2.775—4 

total 

4.675—4 

—3.075—4 

3.675-3 

ratioproj 

N/A 

.25 

4.5 

ratiOquad 

N/A 

1.1 

1.1 

Fig.  5.  Finite  element  solution  components  for  Case  1. 
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Fig.  6.  Adjoint  solution  components  for  Case  1.  Shown  are  plots  for  an  adjoint  solution  using  a  40  x  40  grid  next  to  64  x  64  grid. 
A  solution  on  a  finer  grid  is  used  to  compute  the  estimates. 


4.4  Test  Case  2 

The  next  test  problem  presents  a  more  difficult  solution  for  which  the  geometric  projection  error  is  by 
far  the  largest  source  of  error.  The  grid  is  40  x  40  next  to  64  x  64  and  the  boundary  conditions  are  g  =  0 
on  both  subdomains.  Fig.  7  shows  profiles  of  the  source  and  diffusivity,  while  Fig.  8  shows  the  adjoint 
data.  1 


/GO 


Fig.  7.  Source  /  (left)  and  diffusivity  a  (right)  profiles  for  Test  Case  2.  The  plots  are  shown  in  one  dimension  since  the  source 
and  diffusivity  have  no  variation  in  the  ^-direction. 


Fig.  8.  Adjoint  data  profiles  for  Case  2.  The  plots  of  i f/Ux,  \f/Uy,  and  i ffp  are  shown  in  one  dimension  because  they  have  no  variation 
in  the  ^-direction. 


(*-br 


1  The  shapes  in  Fig.  7  and  8  are  based  on  a  normalized  gaussian  of  the  form  ae  ^2 —  •  The  parameter  a  is  for  normalization 
and  is  set  to  a  —  .  The  parameter  b  determiness  the  location  of  the  peak  and  is  set  to  zero  to  coincide  with  the  interaface.  The 

parameter  c  determines  the  width  of  the  peak  and  is  set  to  c  —  .2.  In  the  case  of  diffusivity  the  dip  is  produced  using  the  function 


10.3—  10*  (! 


,e-{x-br 

2C2 


)• 
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Because  the  source  is  large  but  the  diffusivity  is  small  along  the  interface,  the  solution  changes 
rapidly  near  this  region.  This  leads  to  relatively  large  errors  near  the  interface  for  the  geometric  finite 
volume  method.  When  the  adjoint  data  is  concentrated  near  the  interface,  the  relative  size  of  these  errors 
is  revealed.  Table  8  lists  the  error  terms.  For  this  particular  example  problem,  and  this  particular  error 
measure,  the  error  due  to  geometric  projection  is  nearly  eighty  times  the  total  error  associated  with  the 
residuals.  Fig.  9  shows  the  solution  components  for  the  finite  element  case.  Fig.  10  shows  the  solution 
components  for  the  geometric  finite  volume  case,  and  Fig.  1 1  shows  the  adjoint  solution. 


Table  8.  Error  terms  for  Case  2.  The  forward  grid  is  40  x  40  next  to  64  x  64.  The  adjoint  grid  is  160  x  160  next  to  256  x  256. 


term 

MFE 

GFV  (linear) 

GFV  (constant) 

i 

(Rur-tf-nfK) 

1.6E—5 

1.9/?— 5 

2.4/?— 5 

2 

(Rur-Qr  ~  nhffi) 

-2.6E-5 

—2.6/?— 5 

-2.5E-5 

3 

(Rputt-nzi) 

—3.11?— 5 

-3.1/?— 5 

—3.1/?— 5 

4 

4.6E-5 

4.6/?— 5 

4.6/?— 5 

5 

(Rz,ph-zl,p")ri 

4.81?— 8 

9.4/?— 6 

2.6/?— 5 

6 

(PR-*Ltfh)-Zk,n-nWL)r, 

N/A 

2.3/?— 3 

2.7/?— 3 

7 

<nV,'  -ft-tfW.j.ZWr, 

N/A 

9.0/?— 4 

8.4/?— 3 

8 

N/A 

1.2/?— 3 

1.2/?— 3 

9 

QEu,  (nW) 

N/A 

—2.4/?— 4 

—2.5/?— 4 

total 

5.1/?— 6 

4.2/?— 3 

1.2/?— 2 

ratioproj 

N/A 

25 

73 

ratioquacj 

N/A 

11 

9.7 

Ph  uhx  uhy 

Fig.  10.  Geometric  finite  volume  solution  components  for  Case  2.  Zooming  in  reveals  that  wj  is  discontinuous  across  the  interface. 
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Fig.  1 1 .  Adjoint  solution  components  for  Case  2.  Shown  are  plots  for  an  adjoint  solution  using  a  40  x  40  grid  next  to  64  x  64 
grid.  A  solution  on  a  finer  grid  is  used  to  compute  the  estimates. 


4.5  Test  Case  3 


In  our  final  example,  we  examine  a  problem  that  places  only  one  cell  in  the  ^-direction  in  one  of  the 
subdomains.  Such  a  grid  is  only  appropriate  if  the  solution  in  that  subdomain  is  essentially  one  dimen¬ 
sional,  and  varies  only  parallel  to  the  interface.  This  situation  arises  in  core-edge  coupling  in  a  tokamak 
fusion  reactor. 

We  construct  a  problem  with  a  solution  that  is  very  nearly  one  dimensional  in  one  subdomain,  and 
contains  variation  in  the  second  dimension  well  away  from  the  interface.  The  pressure  component  of 
the  solution  is 


p(x,y) 


=  COS 


(rty± 2) 

l  8 


+  0.3  sin(TTx) 


1-tanh  (2(1.5-;y))‘ 
2 


(4.2) 


The  grid  is  1  x  32  next  to  32  x  32  and  the  boundary  conditions  are  provided  by  evaluating  the  known 
solution  at  the  outer  domain  boundaries.  The  source  for  the  problem  is  computed  by  substituting  the 
chosen  solution  into  the  PDE.  The  diffusivity  a  is  one  in  both  <T2/  and  Qr.  The  adjoint  data  is  concen¬ 
trated  in  the  finer  subdomain,  and  is  shown  in  Fig.  12. 


Fig.  12.  Adjoint  data  profiles  for  Case  3.  The  plots  of  \j/Ux,  y/u  ,  and  y/p  are  shown  in  one  dimension  because  they  have  no 
variation  in  the  ^-direction,  and  i//£  is  a  one  dimensional  function  defined  on  the  interface. 


Table  9  lists  the  error  terms.  For  this  example  problem,  the  contribution  due  to  geometric  projection 
with  linear  extrapolation  is  approximately  ten  times  the  total  contribution  associated  with  the  residuals, 
despite  the  fact  that  the  solution  is  changing  slowly  near  the  interface.  The  projection  contribution  is 
much  larger  if  constant  extrapolation  is  used.  Fig.  13  shows  the  solution  components  for  the  finite 
element  case.  Fig.  14  shows  the  solution  components  for  the  geometric  finite  volume  case,  and  Fig.  15 
shows  the  adjoint  solution  components. 
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Table  9.  Error  terms  for  Case  3.  The  forward  grid  is  1  x  32  next  to  32  X  32.  The  adjoint  grid  is  128  X  128  next  to  128  x  128. 


term 

MFE 

GFV  (linear) 

GFV  (constant) 

i 

3.9E-9 

6.8E-7 

-1.5E-5 

2 

C RuR-K-nhRK ) 

1.8E— 5 

1.8E-5 

1.8E-5 

3 

—4.0E—6 

-4.0E-6 

-4.0E-6 

4 

~PrQ) 

1.2E-6 

7.2E-6 

7.2E-6 

5 

(R^p»-zhp'')ri 

0 

-3.8E-7 

1.7E-5 

6 

{Pr-^l(Pr)  -E,h,n-  nf  ^i)r, 

N/A 

1.5E-4 

-6.4E-3 

7 

{n-ut-PL^R(pt),Z«p'')r, 

N/A 

-6.8E-5 

3. IE— 3 

8 

^uAn'M) 

N/A 

-2.4E-5 

2.5E-3 

9 

QEUR(n£K) 

N/A 

2.0E-5 

1.9E-5 

total 

2. IE— 5 

1.0E-4 

-8.0E-4 

ratioproj 

N/A 

7.3 

155 

ratioqllad 

N/A 

1.5 

41 

Fig.  13.  Finite  element  solution  components  for  Case  3. 
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Fig.  14.  Geometric  finite  volume  solution  components  for  Case  3  computed  using  linear  extrapolation. 
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Fig.  15.  Adjoint  solution  components  for  Case  3.  These  plots  are  the  adjoint  solution  using  64  x  64  next  to  64  x  64  meshes.  The 
estimates  were  computing  using  a  finer  grid. 


5.  Iterative  solvers  and  coupling  strategies 

In  practice,  iterative  solution  of  the  coupled  system  is  often  employed.  The  specific  choice  of  solution 
method  is  often  constrained  by  certain  computational  logistics,  such  as  the  state  of  existing  codes  and 
data  structures.  We  briefly  discuss  some  aspects  of  iterative  solution.  The  primary  goal  is  to  show 
that  iterative  solution  strategies  applied  to  systems  like  (2.10)  can  also  be  applied  to  systems  like  (2.7) 
without  large  changes  to  the  computational  structure.  We  do  not  discuss  the  convergence  of  iterative 
solvers. 


5.1  Iteration  on  the  primary  variable 

A  common  iterative  technique  for  the  geometric  finite  volume  method  (2.10)  is  to  start  with  an  initial 
guess  (pR,pR)  and  proceed  with  the  iteration 


'  AL 

0 

\  H+1 1 

'  FL  ' 

0 

Cd  ' 

Pi, 

0 

Ar 

.  Pit  '  . 

.  . 

Cn 

0 

.  Pr  . 

i =  0, 1,2, .... 


(5.1) 


This  iteration  requires  only  the  inversion  of  ,4/  and  Ar,  that  is,  only  single  domain  component  solves. 
The  application  of  Co  and  Cn  can  be  viewed  as  the  coupling  strategy,  in  which  information  is  swapped 
between  the  subdomains. 

It  is  possible  to  use  an  iteration  of  this  type  on  the  finite  element  system  (2.7)  as  well.  We  must  first 
reduce  to  a  system  in  p  by  a  preprocessing  procedure.  We  first  eliminate  ul  and  ur,  which  results  in 


-  b]ml  1bl  o  -bJm,  'ey 

PL 

0  btrmr1br  -BtrMr1Cr 

PR 

C[M£1Bl  CtrMrxBr  -(■ CTLMllCL+CTRMRlCR ) 

[  %  \ 

fl+btlm-1dl 

Fr+BtrMr1Dr 

CtlMl1Dl+CtrMr1Dr 


which  we  write  succinctly  as 


Gl 

0 


0 

Gr 

HTr 


-hl 

PL 

Rl 

-Hr 

PR 

= 

Rr 

(Kl  +  Kr) 

l  4  J 

Sl  +  Sr 

(5.2) 


(5.3) 
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We  then  eliminate  |  to  obtain 


GL-HL(KL  +  KR)~lH[ 

~Hl{Kl  +  Kr)~1Htr 

PL 

-hr(kl+kr)~1h[ 

GR-HR(KL+KRylHTR 

.  PR  . 

'  Rl~ 

Rr- 

■  (5.4) 


System  (5.4)  has  the  same  structure  as  (2.10),  so  an  iteration  analogous  to  (5.1)  can  be  applied.  The 
stencil  within  the  diagonal  blocks  of  (5.4)  is  very  close,  but  not  identical,  to  the  stencil  of  a  single 
domain  discretization.  The  difference  occurs  only  in  the  stencil  corresponding  to  cells  touching  the 
interface. 

In  some  cases,  e.g.,  the  use  of  black  box  single  domain  solvers,  it  is  necessary  to  construct  a  system 
in  which  the  diagonal  blocks  correspond  exactly  to  single  domain  discretizations.  If  this  is  the  case,  the 
strategy  of  “discretization  consistent  interface  conditions”  provides  a  partial  solution.  In  this  strategy, 
the  diagonal  blocks  are  single  domain  discretizations,  just  as  in  (2.10).  The  off  diagonal  blocks  are 
populated  by  writing  down  both  the  Dirichlet  and  Neumann  boundary  condition  equations  for  every 
cell  touching  the  interface,  rearranging  those  equations  to  isolate  the  boundary  value  terms  and  setting 
those  terms  equal  to  each  other  across  the  interface.  If  the  cell  ratio  along  the  interface  is  integer, 
such  as  4  next  to  8,  the  resulting  system  is  algebraically  equivalent  to  (5.4).  If  the  cell  ratio  is  not  an 
integer  ratio,  such  as  5  next  to  8,  the  equality  of  boundary  value  terms  across  the  interface  can  only  be 
enforced  approximately,  and  the  resulting  system  is  not  exactly  equivalent  to  (5.4).  While  a  complete 
discussion  of  the  implementation  of  discretization  consistent  interface  conditions  is  beyond  the  scope 
of  this  paper,  it  is  worth  consideration  as  an  alternative  to  the  full  mortar  method  in  cases  where  the 
computational  structure  is  constrained  by  black  box  single  domain  solvers  in  combination  with  iteration 
on  the  primary  variables.  The  concept  of  discretization  consistent  interface  conditions  is  similar  to 
strategies  employed  in  Farhat  et  al.  (1998)  and  Edwards  &  Rogers  (1998).  We  should  remark  that  the 
former  paper  recommended  against  mortar  methods  for  the  fluid-structure  interaction  problem,  due  to 
the  lack  of  theory  on  optimal  convergence  and  a  need  to  invert  a  large  interface  matrix.  However,  for 
the  problem  considered  in  this  paper,  the  mortar  method  does  achieve  optimal  convergence.  Moreover, 
we  presented  several  computational  strategies  that  do  not  require  inversion  of  an  interface  matrix. 


5.2  Iteration  on  interface  variables 

An  alternative  iterative  strategy  (Glowinski  &  Wheeler,  1988)  uses  the  interface  variables  as  the  primary 
variables.  If  we  combine  the  u  and  p  variables  into  the  symbol  yr,  then  system  (2.7)  can  be  written  as 


& 

o 

$ 

’  Wl  " 

'  &L 

0  S^R  ^R 

y>R 
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&R 

1 - 

o 

1 _ 

_  4  . 

0 

We  eliminate  yt  as 

y/i  =  #fi-\&i-tf£),  i  =  L,R, 
which  gives  the  following  system  for  t, : 


(5.6) 
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If  a  Krylov  method  is  applied  to  system  (5.6),  then  only  matrix  vector  products  involving  the  matrix 
on  the  left  are  required.  Since  this  matrix  contains  1  and  .s//(,  1 ,  obtaining  a  matrix  vector  product 
amounts  to  performing  single  domain  component  solves.  Once  |  is  obtained,  y/  is  recovered  as  above. 
In  the  setting  of  geometric  coupling,  we  rewrite  the  geometric  finite  volume  system  as 

'  Al  0  UD  0 
0  Ar  0  UN 
0  Ed  -I  0 
En  0  0-7 

where  ,4/  and  Ar  are  single  domain  finite  volume  systems,  and  the  coupling  strategy  by  which  Dirichlet 
(D)  and  Neumann  (N)  data  is  provided  by  the  opposite  subdomain  is  defined  by 

E  \  /)  =  N  and  EqPr  =  D. 

Eliminating  D  and  N  from  system  (5.7)  gives 


Al 

uded  ' 

PL 

'  fl  ' 

UnEn 

Ar 

.  PR  . 

.  F R  . 

which  is  identical  to  (2.10).  If  instead  we  eliminate  pi  and  pr,  the  system  (5.7)  becomes 

7  EDA~lUN 
ENAllUD  I 

which  allows  for  an  iteration  of  the  form  of  (5.1)  on  the  values  D  and  N,  from  which  the  primary 
variables  can  be  recovered.  Solving  (5.8)  by  iteration  is  analogous  to  solving  (5.6)  by  iteration,  and  both 
require  only  component  solves. 

6.  Conclusion 

We  compared  the  accuracy  and  performance  of  two  numerical  approaches  to  solving  systems  of  partial 
differential  equations.  The  equations  were  posed  on  adjoining  domains  which  share  a  common  boundary 
interface  on  which  are  imposed  boundary  conditions.  We  treated  the  important  case  of  different  and 
non-matching  meshes  being  used  on  the  two  domains.  The  first  widely  used  approach  was  based  on  a 
finite  volume  method  employing  ad  hoc  projections  on  the  interface  to  relate  approximations  on  the  two 
domains.  The  second  approach  used  the  mathematically-founded  mortar  mixed  finite  element  method. 
To  quantify  the  performance,  we  used  a  goal-oriented  a-posteriori  error  estimate  that  quantifies  various 
aspects  of  discretization  error  to  the  overall  error.  The  performance  difference  that  we  found  may  not 
be  surprising  in  some  cases.  However,  we  believe  that  there  is  a  perception  in  part  of  the  scientific 
community  concerned  with  multiphysics  systems  that  if  the  solution  is  smooth  near  the  interface,  then  it 
is  not  very  important  exactly  how  the  coupling  is  accomplished.  We  found  that,  on  the  contrary,  that  the 
error  associated  with  ad  hoc  coupling  approaches  may  be  large  in  practical  situations.  The  deterioration 
in  accuracy  was  shown  to  be  due  mainly  to  incorrect  transfer  of  information  (or  projection  error)  across 
the  interface.  Moreover,  we  also  showed  that  mortar  methods  can  be  used  with  black  box  component 
solves,  thus  permitting  an  efficient  and  practical  implementation  of  the  mortar  coupling  approach  within 
legacy  codes. 
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